/* ------------------------------------------------------------------------ */
/* Header File, contains global declarations                                */
/* ------------------------------------------------------------------------ */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define GAMMA           1.4
#define GAS_CONST	    287.0
#define EPSILON         1.0e-15
#define DDK             0.1
#define MINNON          3

/* Some useful functions */
#define max(x,y)        ( ((x) >= (y)) ? (x) : (y) )
#define min(x,y)        ( ((x) <= (y)) ? (x) : (y) )
#define printxy(i)      printf("x = %f,  y = %f\n",p[i].x,p[i].y)
#define sign(x)         ( ((x) >= 0.0) ? 1.0 : -1.0 )
#define SWAP(a,b)		{ temp=(a); (a)=(b); (b)=temp; }
#define LIMIT1(a,b)		( max(0.0, (a*a + a*b + EPSILON)/(a*a + b*b + EPSILON)) )
#define LIMIT(a,b)		( max( 0.0, (2.0*a*b + EPSILON)/(a*a + b*b + EPSILON) ) )
#define LIMIT2(a,b)      ( ((a)*(b)>0.0) ? MYLIM(a,b) : 0.0) 
#define MYLIM(a,b)      ( 2.0*sqrt((a)*(b))/(fabs(a) + fabs(b)) )


#define INTERNAL_FLOW   0
#define EXTERNAL_FLOW   1

#define FORWARD         0
#define BACKWARD        1

#define SCRATCH         0
#define RESTART         1

/* Point types */
#define INTERIOR        0
#define SOLID_WALL      1
#define OUTER_BOUND     2
#define FIRST_LAYER     3
#define KUTTA           4
#define CONST           5
#define SUP_OUT         6
#define INTERPOLATE     7
#define CARTESIAN       8

/* Time stepping */
#define GLOBAL          0
#define LOCAL           1

#define NO              0
#define YES             1

#define FAILURE         0
#define SUCCESS         1

#define QLIMIT          0
#define QDD             1

/* Wall bc */
#define SWBC            0
#define PSEUDO          1
#define NWBC            2

#define WEIGHTe(s,s1)   ( exp(-0.5*s/s1) )
#define WEIGHTm(s,s1)   ( (s<s1)?1.0:(s1/s) )
#define WEIGHT4(s,s1)   (s1/s)*(s1/s)
#define WEIGHT(s,s1)    (1.0/s)
#define WEIGHT3(s,s1)   (1.0/sqrt(s))
#define WEIGHT0(s,s1)   (1.0)

double cfl,wall_cfl,swbc_relax,mach_inf,alpha_inf,beta_inf,temp_inf,rho_inf;
double pre_inf,u_inf,v_inf,w_inf,q_inf,H_inf,S_inf,qinf[4];
int no_grid_points,ntri,rk_stage,NOIMP,LAST_ITER,MAX_ITER,no_spl_points,spl_point[50];
int nsolid,solidEdge[1000][2],**tri;
int GridType,startmode,time_mode,wallmode,qmode,FlowType;
int ResFreq,ResDisp;
int IMAX,JMAX; /* Required if using Cartesian or structured grid */
double MIN_RESIDUE,dt_cfl,time_elapsed;
double Cl,Cd,KKK,MaxCond,MaxCondWeighted,MinEdge,MaxEdge,DsVenkat;
double AAA[3],BBB[3];
double **var_old;
char *InputFile,GridFile[50],ResFile[50],ParFile[50];
char OutFile[50],PointFile[50],VigieFile[50],CpFile[50];

struct node{
	double x,y,z,*t1, *t2, *n;
	double var[5],var_old[5],qx[5],qy[5],qz[5],phi[5];
	double *X, *Y, *Z, dt, dsm, ds2;
	int point_type,total_conn,*Conn;
};

void Line();
double ERRF(double);
int Gaussj(double**,int,double*);
void read_input();
void parameters_used();
void geometric(struct node*);
void output(struct node*);
//void rotationmat(double, double, double, double[][]);
void rotationmat(double, double, double, double**);
//void rotate(double*,double[][]);
void rotate(double*,double**);
//void rotate_back(double*,double[][]);
void rotate_back(double*,double**);
void rotate_x(double*, double*, double*, struct node*);
void lift_drag(struct node*);
void minmax(int,int,double*,double**,double*,double*,double*);
void venkat(int,int,double*,double**,double*,double*,double*,double);
void VanAlbada(int,int,double*,double**,double*,double*,double*);
void interpolate(int point, struct node* p);

void xyzflux(double*,double*,double*,double*);
void xflux_plus(double*,double*);
void xflux_minus(double*,double*);
void yflux_plus(double*,double*);
void yflux_minus(double*,double*);
void yfluxN_plus(double*,double*);
void yfluxN_minus(double*,double*);
void zfluxN_plus(double*,double*);
void zfluxN_minus(double*,double*);

void FluxKinetic(double*, double*, double*);

void QDerv(struct node*,double**);
/* -------------------------------------------------------------------- */
